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Motivated by high interest in the close relation between string theory and black hole solutions, in 
this paper, we take into account the Einstein-Gauss-Bonnet Lagrangian in the context of massive 
gravity. We examine the possibility of black hole in this regard, and discuss the types of horizons. 

Next, we calculate conserved and thermodynamic quantities and check the validity of the first 
law of thermodynamics. In addition, we investigate the stability of these black holes in context 
of canonical ensemble. We show that number, type and place of phase transition points may 
be significantly affected by different parameters. Next, by considering cosmological constant as 
thermodynamical pressure, we will extend phase space and calculate critical values. Then, we 
construct thermodynamical spacetime by considering mass as thermodynamical potential. We study 
geometrical thermodynamics of these black holes in context of heat capacity and extended phase 
space. We show that studying heat capacity, geometrical thermodynamics and critical behavior in 
extended phase space lead to consistent results. Finally, we will employ a new method for obtaining 
critical values and show that the results of this method are consistent with those of other methods. 


I. INTRODUCTION 

Cosmological observations show that about 95% of the universe is made of dark energy and dark matter. Although 
the basic nature of the mentioned dark sector of the universe is an open question, in order to interpret their effects, 
many authors have proposed modified gravity models such as; scalar-tensor theories [l] (or their conformally related; 
F(R) gravity theories @-01), brane world cosmology 0] and Lovelock gravity [!]. This is one of the main motivations 
for considering Lovelock theory. 

On the other hand, regarding higher dimensional gravity, the reasonable effects of higher curvature terms in higher 
derivative gravity theories have been investigated. One of the well-known theories of higher derivative gravity is 
Lovelock theory which is a natural generalization of Einstein theory in higher dimensions. Taking into account the 
first additional term of Einstein gravity in the context of Lovelock theory (Gauss-Bonnet (GB) gravity), it is believed 
that GB gravity can solve some of the shortcomings of Einstein theory Q. GB gravity consists curvature-squared terms 
which, interestingly, is free of ghosts and the corresponding field equations contain no more than second derivatives 
of the metric [§|, H . Another interesting aspect of GB theory is that it can be arisen from the low-energy limit of 
heterotic string theory ESQ- It was shown that the low-energy expansion of a heterotic (closed) string theory 
effective action has GB term as well as a scalar field a The functional form of scalar field remains a fascinating 
and challenging topic of observational research and some authors have stipulated a trivial constant value for the 
scalar field, and therefore, ignore its effect on the solutions. Considering GB gravity context, black hole solutions and 
brane-world model and their interesting behavior have been investigated in literature IT3 | . Recently, it was shown that 
conserved charges of Einstein-Gauss-Bonnet AdS gravity are presented in electric part of the Weyl tensor which gives 
a generalization to conformal mass definition [l3|. Holographic thermalization in GB gravity through the Wilson loop 
and the holographic entanglement entropy has been investigated in Ref. [l4[ . In addition, regarding p -wave phase 
transition, it was shown that although magnetic field has a positive contribution to superconductor phase transition, 
the effect of GB parameter is negative in case of this phase transition [lj|. In the context of AdS/CFT correspondence, 
there has been several studies with the subject of GB gravity [16] . In Ref. iia, Stuckelberg scalar field has been 
employed to study phase transitions in a dual quantum field theory by considering pure GB gravity in asymptotically 
AdS space. On the other hand, in Ref. [18], it was shown that restriction due to causality violation in GB gravity 
points to the need to complete the theory with an infinite tower of massive higher-spin states, as it happens in string 
theory . Motivated by the recent results mentioned above, we study thermodynamic behavior of black holes in 
GB-massive gravity. 

Before we proceed, let us provide a motivation for considering massive gravity. One of the interesting open question 
of gravitational field is graviton and its properties. Although graviton is massless in the Einstein gravity, a natural 
question is whether one can build a self-consistent gravity theory if graviton is massive. Therefore, the Einstein gravity 
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may be modified to satisfy this point, where massive gravity is one such modification. From the theoretical perspec¬ 
tive, the shear difficulty of constructing a consistent theory of massive gravity makes the subject more interesting. 
Therefore, there are limited considerable works in the massive gravity context. The development of the ghost-free 
theory with massive gravitons which are non-interactive in flat background was done in Ref. [2^] ■ It was shown that 
generalization of this theory to curved background leads to existence of the ghost instabilities 21j. The effects of 


quantum interactions of massive gravity and a nonlinear theory of massive gravity in absence of ghost field 1 2 21 1 were 
investigated in Refs. [23!,[24j]. In addition, a special class of charged massive black holes has been investigated 25 1 
(see |26( for more details regarding considering massive gravity). A new class of nontrivial massive black holes in 
AdS spacetime was investigated in [27l [28| . In this theory of the massive gravity, graviton has similar behavior as 
lattice in holographic conductor model. In other words, due to exhibiting a Drude peak which approaches to delta 
function in limit of massless gravity, graviton in this theory plays the role of lattice. Thermodynamical behavior, 
P — V criticality and geometrothermodynamics (GTD) of the Vegh’s massive gravity has been, recently, investigated 
[29| - l3T1 |. Some holographic consequences of graviton mass have been investigated in Ref. 0]. In this paper, we intend 
to generalize Einstein gravity action by adding GB Lagrangian and obtain exact charged solutions of GB-Massive 
gravity. In general, higher derivative gravity terms (such as GB) are corrections that become relevant at high-energy 
regime while the mass terms (of massive gravity) are relevant at low-energy. Here, we do not regard the energy 
regime, directly, and we look for the weight of GB and mass terms contributions into the thermodynamical behavior 
of a typical black hole. 

Obtaining the solutions of GB-massive gravity, we should check their stability. In general, one may categorize the 
stability criteria into two classes; dynamical stability and thermodynamical one. In the present paper, we relax the 
dynamical stability and instead, we focus on thermodynamic behavior. Thermodynamical aspects of the black holes 
are among the popular subjects for studying black hole’s behavior. Thermal stability of the black hole is of great 
importance and has been widely investigated in literature (33j. This is due to the fact that instability of black holes 
puts restriction (conditions) for validation of the physical solutions. One of the approaches for studying thermal 
stability is through canonical ensemble with investigation of the heat capacity. On the other hand, one can look for 
the phase transition points of solutions by calculating root(s) and divergence point(s) of the heat capacity. Also, 
the thermal stability and its corresponding conditions can be employed in studies conducted in context of AdS/ CFT 
correspondence. 

Recently, there has been a renewed interpretation regarding cosmological constant as a thermodynamical variable. 
It was pointed out that considering cosmological constant as a thermodynamical variable may lead to removing 
ensemble dependency in studying stability in context of canonical (heat capacity) and grand canonical (Hessian 
matrix) ensembles for BTZ black holes [34|. In addition, in context of AdS/CFT correspondence, it was shown 
that variation of cosmological constant in black holes corresponds to variation of the number of colors in Yang-Mills 
theory residing in the boundary of the spacetime (35|. There has been a consideration of the cosmological constant 
as a state-dependant parameter in two dimensional clilaton gravity [36[. In an interesting paper, Hawking and Page 
showed an interpretation of confinement/deconfinement phase transition in the dual strongly coupled gauge theory 
for phase transition between the stable large black hole and thermal gas in AdS space [ 37 J. On the other hand, it was 
shown that there is a similarity between phase transition in black holes and liquid/gas phase transition in a Van der 
Waals system [3£|. This consideration has been investigated in literature for different types of black holes [3!|. 

Another method for studying critical behavior of the system is by constructing a spacetime by considering a 
thermodynamical potential and a set of extensive parameters as components of that thermodynamical spacetime. 
The calculated Ricci scalar of constructed thermodynamical spacetime should diverges in place of phase transition 
points that were mentioned in studying heat capacity. There are several approaches for this method that to name 
a few, one can say: Weinhold in which mass is considered as thermodynamical potential with other parameters 
(entropy, electric charge and etc.) as extensive parameters (40], Ruppeiner approach that is built up based on entropy 
as thermodynamical potential | 4l| , it should be pointed out that it was shown that these two metrics are conformably 
related to each other by a conformal factor (temperature) [42| . The other method is the one that Quevedo introduced 
which has mass as thermodynamical potential like Weinhold but with different structure [43|. There are several types 
for Quevedo metric. Recently, it was shown that mentioned approaches in case of some black holes, will not yield 
consistin g re sults with studying heat capacity. In order to overcome this problem, a new metric was introduced in 
Ref. [ 31 L l44l] . The denominator of the Ricci scalar of HPEM metric [E], [44j only contains expressions of numerator 
and denominator of the heat capacity. 

The outline of the paper will be as follows. In Sec. El we introduce action regarding charged massive gravity 
in presence of GB gravity and obtain charged black hole solutions in this gravity. We calculate conserved and 
thermodynamic quantities related to obtained solutions and check the validation of the first law of thermodynamics 
in Sec. ED In Sec. m we study thermal stability of the massive charged GB black hole solutions in canonical 
ensemble using heat capacity. In the next section, considering cosmological constant as thermodynamical pressure, 
we study phase transition of black holes in context of P — V criticality and phase diagrams. In Sec. ED we study 
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thermodynamical behavior of these solutions through geometrical thermodynamics. Sec. I VIII is devoted to calculation 
of critical values in extended phase space by using heat capacity. Last section will be closing remarks. 


II. CHARGED BLACK HOLE SOLUTIONS IN GB-MASSIVE GRAVITY 


The d-dimensional action of GB-massive gravity with a cosmological constant (A) can be written as 


Z = 


1 

167T 


d d xyj~^g 


4 

1Z — 2A + aLcB — T + m 2 ^ CiUj(g, d 1 ) , 

i 


(i) 


where m and TZ are, respectively, the massive parameter and the scalar curvature and F = F^F^ 11 is the Maxwell 
invariant. In addition, a and Lqb are, respectively, the coefficient and the Lagrangian of GB gravity, and dMs a 
fixed symmetric tensor. In Eq. Ci are constants and Hi are symmetric polynomials of the eigenvalues of the d x d 
matrix K= \/g^ a \\> a v 

U x = [AC], 

U 2 = [1C] 2 - [K 2 ] , 

U 3 = [1C] 3 — 3 [1C] [1C 2 ] +2 [1C 3 ] , 

u 4 = [/C] 4 - 6 [1C 2 ] [1C] 2 + 8 [1C 3 ] [1C] + 3 [K 2 ] 2 - 6 [/C 4 ] , 

and L G b = R^sR^ 6 - 4R flv R flv + R 2 where R^ and R^s are Ricci and Riemann tensors. Variation of the 
action m with respect to the metric tensor and the Faraday tensor F^ v , leads to 


G u 


— 


2F ll \F l , 2 


+ mXpv = 0 , 


( 2 ) 


= 0, 


(3) 


where G^ is the Einstein tensor, and x^v are 

= — — (8R pa R^pva — AR p ° x Rv pa \ — 4RR^ U + SR^xR^ + g^LcB ), (4) 

= -j (Wi g„ v - AC m „) - j (i U 2 g - 2 U X TC^ + 2AC£„) - j(H :i g fl , ~ 3 U 2 1C^ 

+6 HxTC% - 6 K\ v ) - y- 4U A F ilv + 12W 2 AC^ - 2411^%, + 24 /C 4 „). (5) 

We consider the metric of d-dimensional metric with the following form 

ds 2 = —f(r)dt 2 + f~ 1 (r)dr 2 + r 2 hijdxidxj, i, j = 1, 2,3,..., n , (6) 

where hijdxidxj is the line element of a (d— 2)-dimensional space with constant curvature (d — 2)(d — 3)fc and volume 
Vd- 2 - We should note that the constant k indicates that the boundary of t = constant and r = constant can be a 
positive (elliptic), zero (flat) or negative (hyperbolic) constant curvature hypersurface. 

Here, we consider the ansatz of Ref. [29j] for the auxiliary reference metric 

d'/i,' = diag( 0, 0, c 2 hij), 

where c is a positive constant. Using the ansatz ([7]), Hi s may be written as [29] 

d 2 c 

1 

r 

d 2 d 3 c 2 
^2 ’ 

d 2 d 3 d A c 3 
^3 ’ 

d 2 d 3 d A d 3 c^ 
r 4 ’ 


U x 

U 2 

u 3 

U A 


( 7 ) 
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in which we used the notation di = d — i. 

I order to study the effects of massive gravity, it is necessary to consider a reference metric, denoted by 4' M „. In 
a special class, one may treat the reference metric in the context of a bimetric theory to credit it with a dynamical 
spin-2 tensor field (see [45| for more details). Another class of considering reference metric comes from the gauge 
invariance property with using of Stuckelberg scalar fields [46l - l48| . It was shown that the number of degrees of freedom 
propagated by any dRGT-type massive gravity in d-dimensional space-time (d > 3) with N scalar fields (Stuckelberg 
fields) equals to |d(d — 3) + N. Using the Stuckelberg language it was shown that the number of propagators is less 
than the number of fields [49|, [5(| . Vanishing of the determinant of the kinetic (Hessian) matrix of the scalar field 
Lagrangian helps us to reduce the number of propagators. Regarding the determinant of the kinetic matrix ( detS ), 
the so-called singular solutions are related to detS = 0 [5lj| . In general, we expect that all regular solutions (detS ^ 0) 
to be tangential to the surface detS = 0 at some point of time. It means that the singular solutions (detS = 0) are 
related to the families of regular solutions. In other words, we can set the conditions in vicinity of the singular surface 
and apply the infinitesimal evolution in time to obtain regular solutions which are tangential to the singular ones 

mm. 

For example, 4-dimensional massive gravity theories admit, generally, up to six propagating modes where five of 
them are physical and one may be a ghost. Depending on the reference metric, the number of degrees of freedom 
that can be propagated, hence the existence of ghost may vary. The existence of the ghost determines whether the 
theory under consideration is actually stable or not. In this paper, the reference metric has the form of 4/ = (0, 0, chij ) 
where in a simpler case for 4-dinrensional spacetime it will have a form of 4/ = (0, 0,1,1). This is a degenerate metric 
which is equivalent to a theory of massive gravity with Stuckelberg fields in the unitary gauge 27}. The reason for 
such consideration is the fact that under a specific coordinate transformation the spatial mass term breaks in spatial 
dimensions while the general covariance is preserved in radial and temporal coordinates. Since this choice of the 
reference metric leads to preservation of the Hamiltonian constraint, one of the degrees of freedoms is removed. In 
addition, the diffeomorphism is not broken in this specific theory of massive gravity. This leads to elimination of 
another degrees of freedom. Therefore, 4-degrees of freedom exist for the propagation with this reference metric. 
Now, considering that only two Stuckelberg fields exist in the diffeomorphism invariant formulation of this theory, 
two degrees of freedom are absent which leads to absence of Boulware-Deser ghost [13, IHJ . The existence of ghost for 
this theory was investigated in Ref. [22j ■ The mentioned reference metric is a singular one. It was shown that this 
theory with any singular metric is a ghost free theory [54j |. In addition, a case of reduced massive gravity with two 
Stuckelberg fields was investigated in Ref. [52j. The extension of the theory to a negative cosmological constant and 
a Maxwell field was done in Ref. [HU ■ The stability of the black brane solutions and the degrees of freedoms that are 
provided by Maxwell field and massive gravity was investigated in detail and it was pointed out that the stability of 
the solutions and degrees of freedoms depends on the choices of the parameters. Regarding the above statements, we 
do not concern for the existence of regular solutions. Therefore in this paper we relax the dynamical stability criteria 
and focus on the thermodynamical properties as well as critical behavior. Use the electrical gauge potential ansatz, 
Ap = H(r)6°, with Maxwell equation 0. we obtain 


H(r) = 


-q 


d 3 r d3 ’’ 


( 8 ) 


where q is an integration constant which is related to the electric charge. Also, the Maxwell equation implies that the 
electric field in d-dimensions is given by 


Ftr ~ r d 2 ' 


(9) 


Now, we want to obtain the static black hole solutions. For this purpose, one may use any components of Eq. 0 
and obtain metric function f(r). One can use different components (tt and X\X\) of Eq. fl3J) which can be written as 


tt = 


d 2 d 3 f 

‘2y-‘2d,2 

[m 2 cci — /'] r 2db / 2 2A r 2d2 2 q 2 


{d 4 d 5 [m 2 c 4 c 4 + af 2 ] r 2di + d 4 [m 2 c 3 c 3 + 2 aff] r 2d V 2 + [m 2 c 2 c 2 - /] r 2d3 

— r 2d-3 k ( ad 4 d 3 [2/ — k] r~ 2 + 2ad 4 f'r ~ 1 — l)}, 


did- 


2«3 


d 2 d 3 


x 4 x 4 = 


2^2d 3 4 { rf 5 rf 6 [m 2 c 4 c 4 + af 2 ] r 2di + d 5 [m 2 c 3 c 3 + 4 aff] r 2dl ' 2 - ^ T — 


+ [m 2 c 2 c 2 + 2 af + 2 aff" — f] r 2d3 


[m 2 cci - 2/'] r 2ds / 2 2 q 2 


d 4 


d 3 d 4 


—r 2d,3 k (ad 3 d 3 [ 2 / — k] r 2 +4ad 5 f'r 1 +2af" - l)}, 


( 10 ) 


( 11 ) 
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FIG. 1: f(r) versus r for A = —1, q — 1, a = 0.3, c = 1.2, ci = 0.6 and d = 5. 

Left panel for k = 0, ci = —1.8, C 3 = 0.5, C 4 = 2, mo = 2, m = 2.100 (dashed line), m = 1.895 (continues line) and m = 1.700 
(dotted line). 

Right panel for k = 1, C2 = 1, C3 = —0.5, mo = 1.5, m = 1.4, C 4 = —0.21 (dashed line), C 4 = —0.18 (continues line) and 
C 4 = —0.16 (dotted line). 


where the prime and double prime are, respectively, the first and second derivatives with respect to r. We can obtain 
the metric function /(r) using the Eqs. (flOl) and (fill) . yielding 


f(r) 

T 


2ad,3d,4 


1 - 


8ad3<i4 

dlC?2 


A + 


did 2 mo q 2 d± 

2 r dl d^r 2d2 


—m 2 diG?2 


d^diC A C4 d3C 3 C3 C 2 C2 


CCi 


2r 4 


2 r 3 


2 r 2 2 


( 12 ) 


where too is an integration constant which is related to the total mass of the black hole. It is notable that, obtained 
metric function m, satisfy the all components of the Eq. 0- Also, in the absence of massive parameter (to = 0), 
the solution CE2D reduces to 


/ ( r ) = k + 


2adsd^ 



Sad^d^ 

d\d2 


A did 2 mo 
2r d i 


Q 2 di 11 

d 3r 2d 2 ~ j ’ 


(13) 


which describes a d-dimensional asymptotically (A)dS topological black hole with a negative, zero or positive constant 
curvature hypersurface. In order to study the effects of the massive gravity on metric function, we have plotted 
following diagrams (Figs. [I] and [2]). 

It is evident that for specific values of different parameters, the metric function will have the usual behavior that one 
can see in GB-Maxwell gravity: there may be two horizons, one extreme horizon and no horizon (naked singularity) 
(Fig. [Tj). Interestingly, considering specific set of values for metric function parameters will lead to modification in 
metric function behavior and variation from the usual behavior. In other words, the black holes in this configuration 
may have following behavior: two normal and one extreme (outer) horizons, four horizons (Fig. [2]). The variation 
in number of the horizons emphasizes the contribution of the massive part. In other words, due to existence of the 
massive part, the geometrical structures of the black holes are modified and hence, their corresponding phenomenology 
is also altered. It should be pointed out that having three or four horizons is a function of the massive coefficients. 

Now, we are in a position to study the existence of the singularity. It is straightforward to show that all curvature 
invariants of the spacetime such as the Weyl square, Ricci square, Ricci and Kretschmann scalars are only the functions 
of /"(r), and thus it is sufficient to study one of them. It is easy to show that 


f'(r) 


2g?2 d\ 


2^3 


f(r) -k \ 
r 2 J 


R aM sR aMS = f" 2 {r) + 2d 2 


2 


(14) 
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FIG. 2: f(r) versus r for A = —1, q = 1, a = 1.4, m = 2, ci = —0.8, C 2 = —0.5, C 3 = 1, C 4 = —0.8, mo = 0.4, k = 1 and d = 5. 
diagrams for c = 1.200 (dashed line), c = 1.082 (continues line) and c = 1.000 (dotted line). 



FIG. 3: Carter-Penrose diagram for the asymptotically adS black holes when the metric function has four real positive roots 
(n < ?'2 <r3 < r±). 


where 


f'(r) = 


■( 1 - 12 ) 2 r 2 


adzdi d\d2$d 


f"(r) = 


1 - 12 


8 rE 


8adsd4r 


2—2 


2r 2 r 


ad^di d\d2&, d\d |12 3 


(15) 

(16) 


in which 


= m 2 diG?2 


2d 3 (i4C 4 C4 3dzC 3 C3 C 2 C2 CCI 


2 r 4 


2 d 2 r 2 


d\d2m.q 2q 2 d\d2 


2r d 


d 3 r 2d 


-3 : 
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FIG. 4: Carter- Penrose diagram for the asymptotically adS black holes when the metric function has three real positive roots 
(ri < r 2 < r 3 ). n in left panel and r 3 in right panel are extreme roots. 




FIG. 5: Carter-Penrose diagram for the asymptotically adS black holes when the metric function has two extreme roots (left 
panel: ri < n) and three real positive roots (right panel: r 1 < r 2 < r 3 ). T 2 in right panel is an extreme root. 


r = 


dd 2 d2irio 

2 r d ~ i 


— mrdid2 


10 d 3 (i 4 C 4 C 4 6 dsc 3 C 3 3c 2 C 2 cci 


+ 


d 2 r 3 


2q 2 d\d2 (2d — 3) 


d 3 r 


2 d, 


— \ /1 + 


8ad 3 d4 

d\d2 


A 


did2m.o q 2 di 
2 r dl d 3 r 2d2 


T 


Since the metric function diverges at r = 0 and is a smooth regular function for r > 0, its first and second derivatives 
are regular functions for r > 0. Therefore, we expect to have only one physical curvature singularity located at r = 0. 
To confirm it, one can investigate the Kretsclimann scalar for various ranges of r. Straightforward calculations show 










that Ra^rysR 01 ^ 15 oc r 2 ( n ^ for small values of radial coordinate. In other words, one concludes 

lim R a B 7 sR a0 ^ 6 = oo, 

r —>-0 


8 

(17) 


which shows that there is an essential singularity located at the origin. For asymptotical behavior of the metric 
solutions, one can series expand the Kretschmann scalar for large values of radial coordinate which leads to 


lim i? Q , / 3 7 < 5 i?“ /37 ' 5 = 


[ d - idid,2 + (4Ag?_i + di ) ad ^ d ^\ — ( d -\ d2 + d ^ d ^ a ) \ Jd \ d2 [ Sd ^ d^aA + d \ d ^\ 

did2d^d\a 2 


(18) 


in which for small values of GB parameter it will yield 


lim R a0lS R a0lS = 


8d -\ A 2 


d \ d , 


4 

d\d2 


A + O (a). 


(19) 


Eqs. (1181) and fflTJl) confirm that asymptotical behavior of the solutions is (A)dS with an effective cosmological 
constant A e ff = A e //(A, a ). 

Here we give a brief discussion regarding Carter Penrose diagram. It is believed that for investigating the conformal 
structure of the solutions, one may use the conformal compactification method to plot the Carter Penrose or conformal 
diagram (see Figs. [3] []] and 0. According to the Carter-Penrose diagrams one finds that, the singularity is timelike 
such as that of Reissner-Nordstrom black holes. In other words, although massive part of metric function may 
change the horizon structure of black holes, it does not affect the type of singularity and asymptotical behavior of 
the solutions. Drawing the Carter-Penrose diagrams, one can find that the causal structure of the solutions are 
asymptotically well behaved. For completeness, we should note that the curvature singularity comes from the nature 
of metric function and the geometrical nature of spacetime. In order to investigate physical singularity, one should 
study the Riemann tensor. It is a matter of calculation to show that the nonzero components of Riemann tensor for 
topological black holes are 


Rtrtr = \ f"(r) 


Rtxitxi = \rf'(r)Q 


R r 


r /'(O p _ _1 ( b ' f ( r )]' 
2 /(r) U - 2 ^ /( r ) 



RxxXiXixi = r 2 [k - /(r)] 0 


where 


0 = 


n— 1 

I ] sin 2 ( Xi ) 

i=l 

1 

n—1 

J] sinh 2 ( Xi ) 
2=1 


k = 1 
k = 0 
k = -1 


and the first and second derivatives of the metric function are given before. Although we cannot analytically check the 
finiteness properties of the Riemann tensor, numerical calculations show that all singularities of the metric function 
(r = 0 and the roots of f(r)) are coordinate singularity except r = 0. The Kretschmann scalar and also the nature of 
Carter-Penrose diagrams confirm that the only physical singularity is located at r = 0. 


III. THERMODYNAMICS 


In this section, we calculate the conserved and thermodynamics quantities of the static black hole solutions in 
d-dimensional GB-massive context and then examine the first law of thermodynamics. 

The Hawking temperature of the black hole at the outer horizon r .(_, may be obtained through the use of the 
definition of surface gravity. Straightforward calculations lead to 


T 


1 

47t efeAf * 372 [Cf + 2fcd3<ii] 


|^2 - 2r^. ^A r 2 ^ 2 + q 2 ^ + d2m 2 ccir"_^ 1/2 j , 


^2 = d2dsd±d5 [afc + to 2 c 4 C4] r^ 2 + d2d% [/c + m 2 c 2 C2] r 2 ^ 1 + d2d3d^m 2 c 3 C3r 2 ^ 3/2 . 


( 20 ) 
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In order to calculate the electric charge of the black hole, one can use the flux of the electromagnetic field at infinity, 
yielding 


Q = 



( 21 ) 


Next, for the electric potential, 17, we use the following definition [jjJ 

U = A, I™ - | r—>r+ = yM. (22) 

d 3 r + 3 

In order to obtain the entropy of the black holes, due to generalization to Gauss-Bonnet gravity, one can use Wald’s 
formula to calculate it pj|. This leads to 


5 = 




2 d 2 d 3 



(23) 


which shows that area law is violated for GB black holes with non-flat horizons (k ^ 0). 

In order to obtain total mass of the black holes, one can use Hamiltonian approach which results into 


M = 


V d2 d 2 m 0 
107 T 


(24) 


Having conserved and thermodynamic quantities at hand, we are in a position to check the first law of thermody¬ 
namics for our solutions. We obtain the mass as a function of the extensive quantities 5 and Q. One may then regard 
the parameters 5 and Q as a complete set of extensive parameters for the mass M(S,Q ) 


M(S,Q) 


2 d 2 (45)* \^A + dnr 2 Q 2 (4S)^ 


2d 

nd\d 2 d 3 (45)^ 


(25) 


where A is 


A = — 2A (45)* + diem 2 (45)* ( Cl + ) + did 2 d 3 m 2 c 3 (45)* ( c 3 + ) . 

V (4 S)*J V (4 S)^J 

We define the intensive parameters conjugate to 5 and Q. These quantities are the temperature and the electric 
potential 


-(f) 


& 


Q 


( dM\ 

~ (w)s' 


(26) 


The results of Eq. EC>1) coincide with Eqs. EOl) and E^l) and, therefore, we find that these conserved and thermo¬ 
dynamic quantities satisfy the first law of black hole thermodynamics with the following form 


dM = TdS+UdQ. 


(27) 


IV. HEAT CAPACITY AND STABILITY IN CANONICAL ENSEMBLE 

In this section, we study the stability of the solutions. In the context of black hole physics, there are two types of 
stability examination that one can employ, the so-called dynamical and thermodynamical stabilities. In this paper, 
we are only interested in thermodynamical stability of the solutions. We investigate thermal stability in context of 
canonical ensemble by calculating the heat capacity. Stability condition states that in order to a black hole being 
thermally stable, its heat capacity must be positive. There are two cases that may happen for an unstable black 
hole: it may have a phase transition and acquire stable state or the obtained solution is not a physical one (no phase 
transition takes place and the system is always unstable). Another important property of the heat capacity that 
motivates one to study it, is phase transition. It is stated that roots and divergence points of the heat capacity 
represent type one and type two phase transition, respectively. In other words, one can employ the numerator and 
denominator of the heat capacity for finding different phase transition points. 
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One can use following relation for calculating the heat capacity 



Considering Eqs. (EU 1 ) and EU 1 ) . it is a matter of calculation to show that heat capacity will be 


Cq = 


B 

C' 


(28) 


(29) 


where 


B = — 


(Zkdzd/ia + r +) 2 d 2 c [i '^ 2 d^d^d^c^ + r^ 5/3 d^d^c^d 1 + r^ 4/3 d^c^c + r+ dl ci^ 


+ 2 d- 2 q 2 ( 2 kd^d^a + r\y 2 + 2Ag?2 ( 2 kd^d^a + r\)~ r. 

—/cc?2c?3 ( 2kdzdiOt + r^_) 2 r^ 2 kd^d^a + r^ 4/3 


2\2 3d 2 /3 


C = [c4<i4C?5 ( 2 kd$dia + 3 r 2 ) c 3 + 2 r 3 G/4C3C 2 — { 2 r\kd^diac 2 — C27 + )c — 4r+fcd4aci] / ^d2d^cm 2 r 2 y i 
—8 q 2 [4/0^7/2^3^40: + 2 d 5 / 2 r + ] r + + 8A [6/0^3^40 + r 2 ] r+ d_1 
+4/od 2 d3 [2/o 2 d 3 d 2 (i5a 2 + fcr 2 d^d^a + r+] r+ dl . 


As one can see, GB parameter and the topology factor (fc) are coupled with each other. It shows that in case of flat 
horizon (fc = 0), the effect of GB gravity is vanished. In other words, in horizon flat solutions, the heat capacity is 
independent of the GB gravity effects. In order to study the thermodynamical behavior of the system, stability and 
phase transition points, we have plotted Figs, [51- [TUI 

It is evident that in case of spherical horizon ( k = 1), there is a root for temperature, rR, in which for r + < cr, 
temperature is negative. Therefore, in this region, obtained solutions are non-physical, tr is a decreasing function 
of massive parameter (Fig. [U] right panel), an increasing function of electric charge (right panels of Figs. [9]) and 
finally independent of GB parameter (Fig. [7] right panel). It should be pointed out that interestingly, increasing and 
decreasing mentioned parameters, will lead to formation of a maximum and minimum. Considering Eq. (1281) . these 
extrema are places in which heat capacity will have divergencies. In case of the topology, interestingly, for k = —1, 
temperature have a root and one divergence point (Fig. [5] right panel) which is not observed in case of horizon flat 
(Fig. [8] middle panel) and spherical symmetric (Fig. [8] right panel). 

For small values of massive parameter, the heat capacity is an increasing function of horizon radius (Fig. [6] left 
panel). For specific values of the massive parameter, a maximum will be formed which is highly sensitive to variation 
of this parameter. For sufficiently large enough of m, there will be two divergencies for the heat capacity and between 
these two divergencies heat capacity is negative, hence, the system is unstable (Fig. [6] left panel). 

Interestingly, the variation of the GB parameter has opposite effect on heat capacity comparing to variation of 
massive parameter. In other words, for small values of the GB parameter, black holes enjoy two divergencies with 
being unstable between these divergencies (Fig. [7] left panel). The interval between these divergencies is a decreasing 
function of GB parameter. For sufficiently small values of a, the heat capacity is only an increasing function of the 
horizon radius without divergency (Fig. [7] left panel). 

As for topological effects, in case of spherical symmetric, no divergency is observed and for the r+ > tr, system 
is in physical stable state (Fig. [5] left panel). On the other hand, in case of horizon flat, two phase transitions are 
observed for heat capacity with black holes being unstable between two divergencies (Fig. [8] left and middle panels). 
Finally, for k = —1, two roots and one divergence point are found. Between these roots and after divergence point, 
system is physical and stable whereas between larger root and divergence point, black holes are thermally unstable 
(Fig. [8] left and middle panels). 

Similar behavior to variation of GB parameter is observed for variation of the charge except for one difference. For 
large values of GB parameter, heat capacity is only an increasing function of the horizon radius whereas by increasing 
the electric charge, the heat capacity will have maximum again (Fig. [H] left panel). 

Finally, in case of dimensions, it is observed that higher dimensionality has contribution to number and place of the 
divergencies for heat capacity. For 5-dimensions, only one root is observed, whereas in case of 6 and 7-dimensions, 
heat capacity has one root and two divergencies which are increasing functions of dimensions (Fig. 11 OH . 

Now we are in a position to discuss phase transitions of the solutions. As it was stated before, roots and divergencies 
of the heat capacity are representing phase transition points of the system. In addition, it should be pointed out that 
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FIG. 6 : Cq (left panel) and T (right panel) versus r+ for q = 1, A = —1 c = ci = C 2 = C 3 = 2, C 4 = 0, a = 0.5, d = 5 and 
fc = 1; m = 0 (continues line), m = 0.35 (dotted line), m = 0.40 (dashed line) and m = 0.50 (dashed-dotted line), ’’different 
scales” 




FIG. 7: Cq (left panel) and T (right panel) versus r+ for q = 1, A = —1, c = ci = C 2 = C 3 = 2, C 4 = 0, m = 0.40, d = 5 and 
k = 1; a = 0.35 (continues line), a = 0.45 (dotted line), a = 0.55 (dashed line) and a = 1 (dashed-dotted line), ’’different 
scales” 


thermodynamical principles state that systems in unstable states go under phase transition to stabilize. Therefore, 
for systems with one root and two divergence points, the following phase transitions happen: for the root there is a 
phase transition from unstable non-physical system to stable and physical one, for smaller divergence point a phase 
transition from larger black holes to smaller ones takes place and finally for larger divergency, systems go under phase 
transition from smaller black holes to larger ones. 


V. P - V CRITICALITY OF CHARGED BLACK HOLES IN GB-MASSIVE GRAVITY 


In this section, we study the phase transition points of the charged black holes in GB-Massive gravity through use 
of P — V criticality and related phase diagrams in spherical symmetric spacetime (k = 1). To do so, we consider 
following relationship between thermodynamical pressure and cosmological constant 


P = - 


A 

8 - 7 T 


( 30 ) 


From thermodynamical point of view, one can indicate that conjugating thermodynamical variable corresponding 
to pressure would be thermodynamical volume. Therefore, in order to calculate the thermodynamical volume of the 
















FIG. 8 : Cq (left and middle panels) and T (right panel) versus r+ for q = 1, A = —1, c = ci = C 2 = C 3 = 2, C 4 = 0, m = 0.40 
d = 5 and a = 0.5; k = 1 (continues line), k = 0 (dotted line) and k = —1 (dashed line), ’’different scales” 




FIG. 9: Cq (left panel) and T (right panel) versus r+ for A = —1, c = ci = C 2 = C 3 = 2, C 4 = 0, a = 0.5, m = 0.40, d = 5 anc 
k — 1; q = 0 (continues line), q = 0.15 (dotted line), q = 0.60 (dashed line) and q = 1 (dashed-dotted line), ’’different scales” 



FIG. 10: Cq versus r+ for q = 1, A = —1, c = ci = ci = C 3 = C 4 = 2, a = 0.5, m = 0.40 and k = 1; d = 5 (continues line) 
d — 6 (dotted line) and d = 7 (dashed line), ’’different scales” 





























solutions, one should use 


13 


V = 



(31) 


Considering cosmological constant as thermodynamical pressure leads to interpretation of mass not only as internal 
energy but also as Enthalpy of thermodynamical system. This interpretation leads to following relation for the Gibbs 
free energy of the system 


G = H -TS = M -TS, 

where by considering Eqs. m, m and (l24l) 

Q _ _ ^0 _ 

lQr d - x d\d 2 d 3 di (r\ + 2a') ’ 


(32) 


(33) 


in which 


2(^2 


0 O = l&'Kdzr 2 ^ (6 d 2 ol + d±r\) P - 4di q 2 [2d 2 d 7 / 2 + ^ 4 ^ 5 / 2 ] - d\d 2 d 3 [2 d 2 a' 2 + d s a'r\ + c^r 2 ] r 

—did 2 d 3 m 2 c \ L ( y d 2 d‘zd±Cid i — d 2 c 2 cr+ — 2cir+) 2 a' + (3d 3 d4C4C 2 + 2d 3 c 3 cr + + c 2 r)}_) djcr 2 ] r 2 ^ 2 , 


and a' = d^d^a. 

By regarding this interpretation and using Eqs. 


, (l30l) and ([3T|) . one can obtain volume as 


y — UJd2 r d i 
d 1 + ’ 


(34) 


which is in agreement with the volume of the spacetime with spherical boundary and radius r + . Because of the 
relation between volume and radius of the black hole, we use horizon radius (specific volume) for obtaining critical 
values. To calculate critical values, one can employ the inflection point properties 


dP 

dr + 


d 2 P 


dr\ J 


1 = 


0 . 


(35) 


As for thermodynamical pressure, by using Eqs. KU1) and (1301) . one can obtain the following equation of state 


P = 


d 2 (2a' + r%) T 


d 2 ( d 5 a' + d 3 r\) m 2 [d 4 d 5 C 4 C 3 + d 4 c 3 c 2 r + c 2 cr^ + Cir 3 ] d 2 d 3 c 


4r 3 


87rr; 


.2 d 2 


167rr 4 


167rr:ji 


(36) 


Now, we are in a position to calculate critical horizon radius. To do so, one should employ Eqs. m and m 
which lead to following relation 


4 d 2 q 2 \pd 7 / 2 a + ( 4 / 2 r+] r 4 = d 2 [l2d 3 a' 2 — I2ar\ + d 3 r 4 ] r 2 + 2 + m 2 [ 6 d 2 d 3 d 4 d 5 C 4 C 4 [2ol + r+) r 


2 d 2 


+3d 2 d 3 d4C 3 c 3 r 2 ^ 1/2 — d 2 d 3 ( 6 a' — r +) 4 c 2 c 2 r+ dl — &d 2 c\ca' r^ x,2 \. 


2d, 


(37) 


It is evident that obtaining critical horizon is not possible, analytically. Therefore, we use numerical method for 
obtaining critical horizon radius and corresponding critical temperature and pressure. The results of this numerical 
evaluation are presented in tables 1 and 2 . 


m r c 

T c 

p Pr.Tc-. 

1 c t 

J c 

0.00000 2.13039 

0.08558 

0.01089 0.27118 

0.10000 1.98493 

0.10527 

0.01415 0.26697 

1.00000 0.91380 

2.63516 

0.77052 0.26719 

2.00000 0.72963 10.78040 3.78491 0.25616 

3.00000 0.65579 24.51754 9.25479 0.24754 


Table (1): q = 1, a' = 0.5, ci = c 2 = c 3 = 2 and d = 5. 
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FIG. 11: P — r+ (left), T — r+ (middle) and G — T (right) diagrams for a' = 0.5, q = 1, m = 0.1, c = ci = ci = cz = 2, a = 0 
and d = 5. 

P — r+ diagram, from up to bottom T = 1.17k, T = T c and T = 0.97k, respectively. 

T — r+ diagram, from up to bottom P = 1.1P C , P = P c and P = 0.9 P c , respectively. 

G — T diagram for P = 0.5P C (continuous line), P = P c (dotted line) and P = 1.5P C (dashed line). 




FIG. 12: P — r+ (left), T — r+ (middle) and G — T (right) diagrams for a' = 0.1, q = 1, m — 0.5, c = ci = C 2 = C 3 = 2, C 4 = 0 
and d = 5. 

P — r+ diagram, from up to bottom T = 1.1T C , T = T C and T = 0.971, respectively. 

T — r+ diagram, from up to bottom P = 1.1P C , P = P c and P = 0.9P C , respectively. 

G — T diagram for P = 0.5P C (continuous line), P = P c (dotted line) and P = 1.5P C (dashed line). 


a' r c T c P c Zp* 

0.00000 0.70743 3.44914 1.82560 0.37444 
0.10000 0.85114 1.70944 0.68803 0.34257 
0.50000 1.19717 0.67256 0.15085 0.26853 
1.00000 1.62237 0.41635 0.05625 0.21922 
2.00000 2.67012 0.26191 0.01729 0.17628 
Table (2): q = 1, m = 0.5, ci = C 2 = C 3 = 2 and d = 5. 


Considering obtained critical values for horizon radius, temperature and pressure, one can plot their corresponding 
phase diagrams (P — r+, T — r + and G — T diagrams). For the economical reasons, we plot phase diagrams for two 
set of critical values (Figs. iTTI and IT2l) . In order to elaborate the effect of variation of GB and massive parameters on 
critical behavior of the system, we plot two sets of comparing diagrams (Figs. fl3l and fT4l). 

The appearance of characteristic swallow tail in G — T diagrams shows that obtained values are critical ones in 
which phase transition takes place (Figs, fill - ITT! right panels). Also, in P — r+ diagrams the existence of region of 
phase transition and critical behavior of the system for different factors of the critical temperature are representing 
phase transition taking place (Figs. |TT| - [□] left panels). As for T — r +1 the formation of subcritical isobars which 
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FIG. 13: P — r+ (left), T — r+ (middle) and G — T (right) diagrams for a' = 0.5, q = 1, c = ci = C 2 = C 3 = 2, C 4 = 0, d » 5, 
m = 0 (continuous line), m = 0.1 (dotted line) and m = 0.3 (dashed line). 

P — r + diagram for T = T c , T — r + diagram for P = P c and G — T diagram for P = 0.5 P c . 



FIG. 14: P — r+ (left), T — r+ (middle) and G — T (right) diagrams for m = 0.5, q = 1, c = ci » C 2 = C 3 = 2, a = 0, d = 5, 
a' = 0 (continuous line), a' = 0.1 (dotted line) and a' = 0.5 (dashed line). 

P — r+ diagram for T = T c , T — r+ diagram for P = P c and G — T diagram for P = 0.5 P c . 


divide the phase states to three different cases (before, in and after phase transition), indicates that calculated critical 
pressure is indeed the one that system goes under phase transition in (Figs. |TT] - [14] middle panels). 

It is evident that, critical horizon radius and universal ratio of are decreasing functions of massive parameter 
whereas critical temperature and pressure are increasing functions of it (see table 1 for more details). Interestingly, 
critical temperature and pressure are highly sensitive to variation of massive parameter. Increasing massive parameter 
for even small value leads to a significant change in the critical temperature. This shows that for black holes to have 
phase transition, it must be heated more significantly. In other words, for acquiring stable state, significant amount 
of energy is needed in order to have phase transition. It is worthwhile to mention that critical horizon also decreases 
highly but not with similar rate as temperature changes. 

Interestingly, except for universal ratio of -Sp 1 , the effects of variation of GB parameter on critical values, are 
completely opposite of massive parameter. In other words, critical temperature and pressure are decreasing functions 
of a whereas critical horizon is an increasing function of it (see table 2 for more details). It is worthwhile to mention 
that the effect of GB parameter is not as significant as the effect of massive parameter. Although both massive 
and GB gravities are extensions for Einstein gravity, their contributions are opposite. The value of GB parameter is 
representing the strength of the gravitational field. Taking this fact into account, one can conclude that strength of 
gravity puts restriction on (weakens) the contribution of the massive gravity. 

Finally, it should be pointed out that subcritical isobars are increasing and decreasing functions of a and m, 
respectively. Subcritical isobars are representing the region in which phase transition takes place (Figs. [13] and [TJ] 
middle panels). Therefore, their increments/reductions will decrease/increase single state region of one phase which 
in our case it is whether small or large black holes. Interestingly, the gap between energy of two different phases is 
an increasing function of both GB and massive parameters. The only differences is that Gibbs energy level of two 
phases are increasing functions of massive parameter (Fig. [13] right panel) and decreasing functions of a (Fig. [14] 
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FIG. 15: 1Z versus r+ for q = 1, A = —1 c = ci = C 2 = C 3 = 2, C 4 = 0, a = 0.5, d = 5 and k = 1; m = 0 (continues line), 
m = 0.35 (dotted line), m = 0.40 (dashed line) and m = 0.50 (dashed-dotted line), ’’different scales” 



FIG. 16: 1Z versus r+ for q = 1, A = —1, c = ci = C 2 = C 3 = 2, ca = 0, m = 0.40, d = 5 and A; = 1; a = 0.35 (continues line), 
a = 0.45 (dotted line), a = 0.55 (dashed line) and a = 1 (dashed-dotted line), ’’different scales” 


right panel). 


VI. GEOMETRICAL THERMODYNAMICS FOR HEAT CAPACITY AND P - V CRITICALITY 


In this section we will study thermodynamical behavior of these black holes through geometrical thermodynamics 
(GTs). As it was pointed out earlier, there are several approaches for constructing thermodynamical spacetime. In 
this paper, we follow HPEM method for studying phase transition points of the heat capacity and later another metric 
for phase transition points of the extended phase space. The HPEM metric has the following form fU [HJ 


ds 


2 

New 


SM s 


f^jrn d 2 M \ 

t i=2 dx'i ) 



( 38 ) 


where Ms = dM/dS , M$s = d 2 M/dS 2 and \i (Xi 7 ^ S) are extensive parameters. For economical reasons, we only 
plot corresponding figures for variation of a and m for studying the heat capacity (Figs. [6] and [7]). By employing Eqs. 
m , (Esd and (l25l) and HPEM metric (Eq. (l38l) h one can plot following diagrams (Figs. fl5l and fl6l). 
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Studying GTs diagrams shows that Ricci scalar of considered metric diverges in places of phase transition points. 
In other words, divergencies of the Ricci scalar coincide with divergencies and roots of the heat capacity (compare 
Figs. H5l and flfil with iGl and 171) . Taking a closer look at the divergencies of the Ricci scalar, one can see that the 
behavior of the Ricci scalar around divergence points can be categorized into three groups: there is a change of sign 
around divergence point in which correspondingly there is a change of sign in heat capacity, divergency toward +oo 
where there is a phase transition from larger black holes to smaller ones and finally divergency toward — oo in which 
system goes under phase transition of smaller to larger black holes. This property of the Ricci scalar enables one to 
recognize the type of phase transition without studying heat capacity. 

In Ref. [55( | a new thermodynamical metric for studying critical behavior in case of P— V criticality was introduced. 
It was also shown that GTs, P — V criticality and heat capacity will lead to consistent results. In other words, number 
and places of phase transitions in these three pictures are uniform and they yield same results. There are two forms 
of GTs metric for studying critical behavior in extended phase space that are given as [ssj 


ds 2 


S (- M ss dS 2 + M QQ dQ 2 + dP 2 ) Case I 
S^(-M ss dS 2 + M QQ dQ 2 + MpdP 2 ) Case II ' 


(39) 


In this paper, we only consider Case I metric for studying GTs in extended phase space. Considering Eqs. (1241) . 
m, USED and ( 1551 ) with obtained critical pressure in table 1 , we plot following diagram (Fig. 1171 ) . 

In plotted graphs, a divergency is observed for Ricci scalar which is due to existence of the root for heat capacity 
(Fig. [IT] left and down-right panels). For pressure smaller than critical pressure, two divergencies for heat capacity 
and Ricci scalar observed (Fig. [FT] up right panel). This behavior is similar to that was observed in Fig. [G] middle 
panel for pressures smaller than critical pressure. In other words, similar to T — r+ diagrams, for pressures smaller 
than critical pressure system will have more than one phase transition point. Next, for pressure being critical pressure 
only one divergence point exists (Fig. [IT] down middle). The place of divergency is exactly where critical horizon 
radius was found. In other words, in case of considering critical pressure, the place of phase transition in extended 
phase space and divergencies of heat capacity and Ricci scalar coincide with each other (compare Fig. [T7]down middle 
with Fig. [6] middle). Finally, for pressures larger than critical pressure no divergence point for heat capacity and 
Ricci scalar was observed which is consistent with result of T — r+ diagram in Fig. [G] 

According to the results, the GTs metric proposed for extended phase space is providing an effective method which 
its results are consistent with thermodynamical concepts. Also, this metric is able to point out phase transition 
points which one can obtain through studying T — r+ and P — r+ diagrams. More importantly, it is evident that heat 
capacity (its divergencies), GTs (its divergencies) and phase diagrams (P — r + , T — r + and G — T) lead to consistent 
results. 


VII. HEAT CAPACITY AND CRITICAL VALUES IN THE EXTENDED PHASE SPACE 


The final section of this paper is devoted to calculation of critical values in extended phase space by using heat 
capacity. In last section, it was seen that critical values in which phase transition takes place in extended phase space 
present themselves as divergencies in heat capacity. Motivated by this result, in Ref. |55| a new method for obtaining 
critical values was introduced which is based on denominator of the heat capacity. 

In this method, one consider denominator of the heat capacity and rewrite it with relation between thermodynamical 
pressure and cosmological constant. Next, one solve the denominator of the heat capacity with respect to pressure. 
Obtained relation for pressure is different from the one that is obtained by using temperature. If obtained relation for 
pressure has maxinrum(s), then in the place of that maximum (horizon radius and pressure), phase transition takes 
place. In other words, instead of mentioned method for obtaining critical values in section V, one can easily study 
the maximums of obtained relation for pressure through denominator of the heat capacity. The advantages of this 
method are pointed out in Ref. [55 ]. 

In order to obtain mentioned relation for pressure, one should use denominator of Eq. (1551) with Eq. (1301) which 
leads to 


in which 


167T r 2 ^ (6a' + 7 +) ’ 


(40) 


(d^d^c^c 3 — d^,C2r\ — 2cir+) 2cfoc/cr^ 2 + (Sd^dsc^c 2 + c^,r- 
-4 q 2 [2 d 7 / 2 a' + d 5 / 2 r+] r+ + cd 2 [2 d 5 a' 2 + dga'r 2 + + d 3 r+] r 2 +' 
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FIG. 17: 1Z (continuous line - Case I), Cq (dashed line) versus r + for q = 1, c = ci = C 2 = C 3 = 2, ca — 0, a — 0.5, d = 5, 

k = 1 and m = 0 . 1 . 

up: P = 0.9P C , ’’for different scales”, 

middle: P = P c ’’for different scales”, 

down: P = 1.1P C , ’’for different scales”. 


Considering mentioned values in tables 1 and 2, we plot following diagrams for Eq. (l40l) . (Fig. fl8l) . 

It is evident from plotted graphs that the maximum of Eq. (l40l) is exactly located at the place in which black hole 
goes under phase transition (compare Fig. (1181) with two tables) with same critical pressure. In other words, the 
maximunrs of Eq. (1401) for different values of parameters, are representing critical horizon radii and pressures in which 
phase transitions take place. Therefore, by employing this method, one is able to calculate critical horizon radius 
and pressure at the same time without going through trouble of obtaining relation for calculating critical horizon 
radius and other relations. Another important property of this method is the consistency of pressure’s behavior with 
thermodynamical concept. In other words, for pressures smaller than critical pressure, two values for pressure of 
phase transition are observed (corresponds to two divergencies for heat capacity in Fig. [17] and T — r + diagram in 
Fig. EU and for pressures larger than critical pressure no pressure for phase transition point is observed. 

In conclusion, one can see that this method enables us to obtain critical values and thermodynamical behavior of 
the system faster and its properties are consistent with thermodynamical concepts. 
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FIG. 18: P versus r + diagrams for q = 1, c = Ci = C 2 = C 3 = 2, C 4 = 0, d = 5 and k = 1. 

left panel: a = 0.5 and m = 0.1 (continues line), P = 0.01415 (dashed line) and r+ = 198493 (dashed-dotted line) 
right panel: a = 0.1 and m = 0.5 (continues line), P = 0.68803 (dashed line) and r+ = 0.85114 (dashed-dotted line). 


VIII. CLOSING REMARKS 

In this paper, we have studied charged massive black holes in Gauss-Bonnet gravity. We obtained metric function 
for this gravity and showed that due to contribution of the massive part, the geometrical structure of black holes 
which includes number of the horizons and their places, have been modified. Therefore, the phenomenology differers 
completely from usual charged Gauss-Bonnet black hole. Plotted Carter -Penrose diagrams showed that although 
massive part of metric function can change the horizon structure of black holes, it does not affect the type of singularity 
and asymptotical behavior of the solutions. Next, we obtained conserved and thermodynamical quantities and showed 
that for obtained values the first law of thermodynamics is valid. 

We also investigated thermal stability of the solutions in context of the canonical ensemble. We showed that 
variation of different parameters, affects the stability conditions of the black holes. Interestingly, it was found that the 
variation of massive and GB parameters has opposite effects on stability and phase transitions of the solutions. Also, 
it was shown that dimensionality has strong contribution to type and number of phase transition which resulted into 
modification of stability conditions. In other words, in specific dimensions, the black hole may only enjoy one type of 
the phase transition whereas in the other dimensions, it may have two types of phase transition in several points. 

Next, by considering cosmological constant as thermodynamical pressure, we study phase transition of these black 
holes in the context of extended phase space. It was shown that thermodynamical volume is independent of GB and 
massive extensions and it only depends on topology of the solutions. Plotted phase diagrams showed that obtained 
critical values are the ones in which phase transition takes place. It was pointed out that interestingly, the effects of 
massive and GB parameters are opposite of each other. The critical temperature and pressure were highly sensitive 
to variation of massive parameter. In opposite, the GB parameter, hence strength of gravitational field put restriction 
on this rapid growth of temperature which will be a controlling factor for massive included systems. 

In addition, through the use of GTs method, phase transition of these black holes was investigated in context 
of heat capacity and extended phase space. It was found that both employed geometrical metrics for studying 
thermodynamical behavior of the system yield consisting results with heat capacity and extended phase space. In 
other words, the divergencies of the Ricci scalar for these metrics coincided with divergencies of the heat capacity and 
thermodynamical behavior of the system in context of phase diagrams. Therefore, these three pictures yield consisting 
results. 

Finally, a new method that was introduced in Ref. (55~| employed for calculating critical pressure and horizon radius. 
It was shown that the results of this method (using maximum of obtained relation for pressure) and study conducted 
in extended phase space were in agreement and obtained critical pressure and horizon radius were the same. 

It will be interesting to study the holographic aspects of obtained solutions in context of superconductors. Also, it 
will be worthwhile to study casual structure and casuality conditions of these systems. In addition, investigation of 
the causal structure of massive gravity in cosmological background by using of interesting Carter- Penrose diagrams 
is attractive. Moreover, following the approach of p52l [S3| . one can study the stability and physical singularities of 
the solutions with massive degrees of freedom. We left these problems for the future works. 
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